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ABSTRACT 

In clusters of galaxies, the specific entropy of intracluster plasma increases 
outwards. Nevertheless, a number of recent studies have shown that the intra- 
cluster medium is subject to buoyancy instabilities due to the effects of cosmic 
rays and anisotropic thermal conduction. In this paper, we present a new numer- 
ical algorithm for simulating such instabilities. This numerical method treats the 
cosmic rays as a fluid, accounts for the diffusion of heat and cosmic rays along 
magnetic field lines, and enforces the condition that the temperature and cosmic- 
ray pressure remain positive. We carry out several tests to ensure the accuracy 
of the code, including the detailed matching of analytic results for the eigenfunc- 
tions and growth rates of linear buoyancy instabilities. This numerical scheme 
will be useful for simulating convection driven by cosmic-ray buoyancy in galaxy 
cluster plasmas and may also be useful for other applications, including fusion 
plasmas, the interstellar medium, and supernovae remnants. 

Subject headings: methods: numerical, conduction, diffusion, convection, MHD, 
plasmas, cosmic rays, instabilities, galaxies: clusters: general, cooling flows 



Introduction 



The hierarchical model of galaxy formation succesfuUy predicts the evolution of baryons 
in the u niverse over a wide range of scales assuming that supernovae feedback is taken into 



account (iKauffmann et al.lll999 
2003 : Springel fc HernquistI 



2003 



Rasera fc Teyssiei 



Somerville fc Primacklll999l : ICole et al.ll2000l : iHatton et al. 



20061 ) . However, the baryon budget 



remains inaccurate in large scale structures (large galaxies, groups or clusters) where the 
total amount of cold gas and stars is overestimated. This overcooling problem is particularly 
critical for galaxy clusters, in which the cooling time near the center of a cluster is often 
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much shorter than a cluster's age. In the absence of heating, one would expect cooling flows 
to form in these clusters, with large amounts of plasma cooling and flowing in towards the 
center. However, the star formation rate in cl uster cores is typ ically 10-100 times lower 
than the predictions of the cooling- flow model (IMcNamaral |2004| ) . and line emission from 
plasma at temperatures lower than one third of the virial temperature of the cluster is weak 
( jPeterson et al.ll2003l : lMcNamarall2004l ). The inconsistency between the cooling-flow model 
and these observations is known as the "cooling-flow problem" . 

A promising hypothesis to solve this puzzle is heating by active galactic nuclei (AGN) 
in cluster cores. Two main arguments support this idea. First, AGN power is expected 
to be a decreasing function of the speciflc entropy at a cluster's center and therefore tends 



naturally towards a se lf-regulated state in which heating balances cooling (INulseru 12004 



Boehringer et al 



2004] . Second, almost all cooling core clusters possess active central ra- 



dio sources Jsilekl I2OO4J ). However, one important problem remains: how is AGN power 



transferred to the ambiant plasma? Over the last decade, a number of n umerical simula 



tions have been carried out to answer this qu e stion. In the flrst simula tions (jChurazov et al. 



2OOII : iQuilis et al.l I2OOII : ISriiggen et al.l I2OO2I : iBriiggen fc Kaiserl I2OO2I ) , thermal energy was 
injected near the center of a 2D or 3D cluster-like hydrostatic proflle. This resulted in hot 
and underdense bubbles, which then rose buoyantly. By agitating the surrounding medium, 
these bubbles were able to reduce the cooling while achieving some correspondence with 
the obse rvations of X-ray c avities seen in roughly one-fourth of the clusters of the Chandra 
archive (IBirzan et al.l 12004 ). Subsequent sir aulations extended these e arlier works to include 
new physical ingredi e nts, such as viscosity (Ruszkowski et al. 2004al lbl: [Reynolds et al. 2005 



Briiggen et all I2OO5I : ISiiacki fc Springell booeh . It was found that viscous dissipation con- 
tributed to the ener gy transfer, and that viscos i ty helped to prevent bubbles from breaking 



up. Other studies (IRevnolds et al. 



Cattaneo fc Teyssiei 



2006 



Heinz et al. 



20011. 12OO2I : lOmma fc BinnevI I2004J : lOmma et all 12004 
20061 ) injected not only thermal energy but also ki 



netic energy in subrelativistic bipolar jets. This approach also leads to cavities, but the 
dynamics are different than in the previous works because of the initial momentum of the 
bubbles and because the energy is deposited over a more narrow range of angles. In this 
context, the importance of turbul ence, magnetohydrody r iamic s effects, and plasma transport 
processes has been underlined by lVernaleo &: Reynoldsl (120061 ) . who suggested that these in- 
gredients could prevent the heating from being highly concentrated along the jet axis, as is 
the case for one-fluid pure-hydrodynamics simulations of jets in clusters that are initially at 
rest. 

The above simulations treated the intracluster medium (ICM) as a single fluid. In 
single-fluid simulations, when AGN-heated plasma at temperature Thot mixes with ambient 
intracluster plasma at temperature Tq, the result is a Maxwellian plasma with a temperature 
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intermediate between Tq and Thot- Although this approach is vahd in clusters if Thot is not 
too large, it breaks down if the hot particles are relativistic or transrelativistic, because 
then Coulomb collisions do not have sufficient time to bring the hot particles into thermal 
equilibrium with the ambient intracluster plasma. If we focus on hot protons, the type of 
collision that brings such protons most rapidly into thermal equilibrium with the background 
plasma is collisions with bac kground electro ns. The time scale for thermal electrons to remove 
energy from a hot proton is iGould I (119721 ) , 



(7 - \)mp'meVp(? 



In 



2meCVpP 



h{AiTe'^ne/me 



2c2 



(1) 



with He the electron density, 7 the Lorentz factor, Vp the proton velocity, e and me the 
electon charge and mass, p and rrip the proton momentum and mass, and h the reduced 
Plank constant. 

For a typical proton energy of ~ 1 GeV (transrelativistic regime) and a typical 
cluster-core electrons density of Ue = 0.01 cm~^, the thermalization time scale is ~ 7 Gyr, 
which is much larger than the time for protons to escape the cluster core via diffusion or 
convection. In this case, the ICM is essentially a two-fluid system similar to the interstellar 
medium of the galaxy, with a thermal background plasma plus a population of high-energy 
particles (cosmic rays). 

There are a few problems with treating a mix of cosmic rays and thermal plasma as a 
single Maxwellian fluid. One is that the single-fluid approximation to the temperature con- 
tains the cosmic-ray contribution to the energy density, and thus overestimates the actual 
temperature of the thermal plasma. If the cosmic-ray energy density is a significant fraction 
of the total energy density, the single-fluid model is unable to accurately predict the temper- 
ature profile of a cluster. In addition, since the thermal conductivity depends sensitively on 
the temperature {kt oc T^^^), and since conduction can m ake an important contribution to 
the heating of a cluster core (jZakamska &: Narayaru |2003| ) errors in the temperature profile 
can also lead to significant secondary errors in the thermal balance of the ICM. 

A more subtle difficulty in applying a one-fluid model to a cosmic-ray /thermal-plasma 
mixture concerns the convective stability of intracluster plasma. It turns out that a radial 
gradient in the cosmic-ray energy density is much more destabilizing than a radial gradient 
in the thermal plasma energy density when the plasma mass density decreases outwards (see 
Eq EH] below). A correct accounting of the fraction of the total pressure contribution by 
cosmic rays is thus essential for understandi ng the convective stability of clusters. A more 
extensive discussion of this point is given by IChandran &: Dennisi ( 120061 ). 



A more accurate treatment of the ICM, in which the cosmic rays are treated as either 
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a second fluid or as collisionless particles, is thus needed. In this paper, we present a new 
numerical algorithm for simulating the ICM that treats the ICM as a two-fluid (cosmic-ray 
plus thermal-plasma) system. We also present the results of a suite of tests for our code. Our 
numerical approach is similar to that of Mathews & Brighenti (2008), who carried out two- 
fluid simulations of cosmic-ray bubbles in the ICM. However, in contrast to this latter study, 
we take thermal conduction and cosmic-ray diffusion to occur almost entirely along mag- 
netic field lines (cross- field transport arising only from numerical diffusion). Such anisotropic 
transport arises in clusters because the Coulomb mean free paths of thermal particles in clus- 
ters are much larger than their gyroradii, and the scattering mean free paths of cosmic rays 
are much larger than their gyroradii. The effects of magnetic fields on conduction are some 
times taken into account in when considering thermal conduction over length scales much 
l arger than the correlation length of the (tangled) intracluster ni agnetic field, — 1 — 10 kpc 

In this case, the conduc- 



dKronbereJ Eooi : IXavlor et al.l boOll. I2OO2 I: IVogt fc E nfflinI 



200 



3.. 



tivity K,T is effectively isotropic ( Rechester fc RosenbluthI 

with a val ue that is ^ 0.1 — 0.2 times the Spitzer thermal condu c tivity for a non-ma. gnetized 



19781 : IChandran fc Cowlevlll998h 



plasma (INaravan fc MedvedevI I2OOII : IChandran fc MaronI l2004l : iMaron et all 120041 ). How- 
ever, on scales < /b, the anisotropy of the therm al conductivity has a powerful effect on the 



convective stability of the intracluster medium jBalbusl[200ol. 



Chandran fc Dennisll2006l : iParrish &: Stonell2007l : iParrish fc Quataert 



2001 



Parrish &: Stone 



20081 : lQuataert1l2008h 



2005 



in such a way as to make convection much more likely than when the conductivity is treated 
as isotropic. This is true even if the magnetic field is so weak that the Lorentz force is 
negligible. In order to simulate buoyancy instabilities and convection in clusters, it is thus 
essential to incorporate anisotropic transport. 

The remainder of this paper is organized as follows. In section [2] we present the basic 
equations of our two-fiuid model. In section [3] we present the total-variation-diminishing 
(TVD) code that we use to solve these equations as well as several numerical tests, focusing 
on the case in which there is no conduction or diffusion. In section H] we present the standard 
numerical discretization method for anisotro pic conduction. We show how it can lead to 
negative temperature as emphasized before by lSharma fc HammettI (120071 ). We then describe 
our new method that does not suffer from negative temperature problems. Tests such as the 
circular conduction test and Sovinec-test are also presented. Finally, in section O we present 
results for the linear buoyancy instabilities involving cosmic rays and anisotropic transport 
and compare our numerical solutions to analytic results. 
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2. Two-fluid equations with anisotropic transport 



In order to carry-out realistic cluster simulations one has to deal with an impressive 
list of components (dark matter, plasma, cosmic rays, magnetic field, stars, supernovae, 
supermassive black holes) and physical ingredients (advection, shocks, induction, gravity, 
anisotropic transport, cooling, energy injection from the AGN, feedback from supernovae, 
jets, viscosity...). In this paper, rather than attempting to simulate all of these processes, we 
focus on developing an accurate and efficient numerical algorithm for simulating collisional 
plasmas pervaded by collisionless cosmic rays. A complete description of the cosmic rays in 
such a system would require us to solve for the cosmic-ray distribution function f{r,p,t), 
where r is the position coordinate, p is momentum, and t is time. The resulting system 
of equations is much more difficult to solve numerically than a system of fluid equations 
because / depends on three momentum coordinates as well as position and time. However, 
in many situations of interest , f is nearly is o tropic in momentum space and can b e treated 
as function of only (r, |p|,t) Jskillinel EoTsl ) . iMiniatil J200lh : iMiniati et all J200l[ ) took ad- 
vantage of this fact with a numerical code, COSMOCR, that solves for the evolution of / 
as a function of both r and p. In this pap er, we adopt th e more simplistic and less com- 
putationally intensive fluid-like approach of iDrury fc Voelkl (Il98ll ). which does not attempt 
to solve for the momentum dependence of /, but instead solves direc tly for the evolu t ion o f 
the cosmic-ray pressure as a function of r an d t. This mode l of iDrurv fc Voelkl (jl98ll ) 



by 


Mathews & 


Briehenti 


Rvu et al. 


(2003 


). In this 



conduction. The resulting equations can be written. 



^ + V.(pt;) = (2) 

dpv „ f BB\ 

+ V. pvv + ptot - — — = pg (3) 



dt ' -y^- ' 4^ 

-^-Vx(«xB) = (4) 

^±^V.[{e\ptot)v-^^^^ = pv.g + V.{K.VT) + V.{D.Ve,r) (5) 

"■+V.(e„.'y) = -p^rV.v + V.{D.Vecr), (6) 



dt 

with the 9 main variables, p the plasma density, pv the plasma momentum, B the magnetic 
field, Ccr the cosmic ray internal energy, and e = 0.5/w^ -|- Cth + ^ct + 0.5-B^/47r the total 
energy. Intermediate variables are eth, the internal thermal energy and ptot = {l — ^)(^th + 
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(7cr — l)ecr + 0.5-B^/47r the total pressure with 7 and 7cr the adiabatic indices of gas and 
cosmic rays, g is an external gravity field (the large scale gravitational potential is mostly 
dominated by stars and dark matter in galaxy clusters). Finally, D and k are the diffusion 
and conduction tensor, which are described further in section m 

In this model, the cosmic rays fiow at the same speed as the thermal plasma, since 
both are frozen to the same magnetic field lines and since wave-particle interactions limit 
the relative motion between cosmic rays and thermal plasma in the direction of the magnetic 
field. On the other hand, because the pitch-angle scattering associated with wave-particle 
interactions is of finite strength, the cosmic rays can diffuse with respect to the thermal 
plasma. We have taken this diffusion, as well as the conduction of heat, to occur entirely 
along magnetic field lines. This is a reasonable approximation in clusters of galaxies, because 
the gyroradii of thermal particles are much shorter than their coUisional mean free paths, 
and the gyroradii of cosmic rays are much shorter than their scattering mean free paths. 
The value of the cosmic-ray diffusion coefficient in clusters of galaxies is not well known. In 
this paper, we use a value D\\ = 10^^ cm^s~^ comparable to the parallel diffusion coefficient 
of 1 GeV protons in the interstellar medium of our galaxy. 

We assume that protons dominate the cosmic-ray energy density in clusters, as is the 
case in the Galaxy. For protons in clusters, the energy loss times associated with Coulomb 
interactions and inelastic collisions (pion production) are typically longer than the growth 
times of the instabilities that we focus on in this paper. Thus, we neglect Coulomb losses 
and pion production in this paper. We note that we also do not include self-gravity, since it 
is not important in the hot intracluster medium. 

It can be seen from equation [6] that the cosmic rays are treated as a fiuid with adiabatic 
index 7^. Thus, if V ■ < at some location, the converging fiow acts to increase the 
cosmic-ray pressure. If we were to model the cosmic-ray distribution function f{r,p,t) as a 
power law in momentum of the form f oc with a l ow eii ergy cutoff, then the effective 



value of 7cr is given by equation [13] of lJubelgas et al.l (120081 ). In this paper we make the 
simple choice that 

7cr = 4/3, (7) 

corresponding to the limit in which a approaches 2 from above, and in which the cosmic-ray 
energy density is dominated by ultra-relativistic particles. 



We note that our approach is in some ways similar to the model of lJubelgas et al.l (120081 ) . 
who incorporated cosmic rays into hydrodynamical simulations of galaxy formation based on 
smoothed particle hydrodynamics. Their approach, like ours, employs an effective adiabatic 
index for the cosmic rays and avoids solvin g for the full cosmic-ray momentum distribution 



function. However, lJubelgas et al.l (120081 ) also develop a framework for incorporating a 
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number of effects that are not treated h ere, including io n izatio n losses, radiation losses, and 
shock acceleration. On the other hand, iJubelgas et al.l (120081 ) assume an isotropic cosmic- 
ray diffusion coefficient, whereas anisotropic transport of both cosmic rays and heat plays a 
central role in our model as well as the buoyancy instabilities that we simulate in section [531 



3. Numerical implementation and tests in the absence of transport 



In this section, we set k, = D = 0. Our numerical method for solving the magnetohydro- 
dynamic-like fMHD-li ke) two-flu i d equ ations is based on the Total Variation Diminishing 
(TVD) MHD code of iPen et al.l (120031 ) which has the advantage of being fast, simple and 
efficient. This TVD MHD code is fully describ ed in 3 papers. Th e appendix of iPenI (119981 ) 
presents the relaxed TVD method that is used. iTrac fc PenI (120031 ) shows the different meth- 
ods for hydrodynamics solver and iPen et al.l (120031 ) describes the MHD code itself. We will 
here recall the main characteristics of this code but the reader should refer to the above 
articles for more complete explanations. 

The fluid solver is a conservative, second-order (in space and time), dimensionally split, 
TVD, upwind scheme. In this relaxing TVD method, each hyperbolic conservation law 
is replaced by a left and a right advection problem with an advection speed called the 
"freezing speed" . By taking this freezing speed equals to the largest eigenvalue of the system 
c = max{\v\ + Cs) (with Cs the sound speed and v the velocity along the updated direction), 
it ensures the scheme to be TVD. The advection problem is then solved using Van-Leer slope 
limiter to reach second order in space and Runge-Kutta integration to reach second order in 
time. 

The magnetic fleld is updated separately in advection-constraints step. A staggered 
grid is used with B deflned on cell surfaces (see Fig{T]) in order to satisfy the divergence-free 
magnetic field condition at machine precision. The advection step is computed using the 
same TVD method as in the fluid solver. This is however easier since the velocity is assumed 
to be flxed (operator splitting). The second order flux is then directly re-used to compute 
the constraint step. Here again Runge-Kutta is used for second order temporal accuracy. 

This method is very efficient because it doesn't need to solve the whole Riemann problem 
and therefore compute all eigenvalues. It only needs the computation of the largest one for 
the freezing speed. The resolution of slow waves is slightly degraded, however the code could 
still resolve shocks using only a few cells. 



The fluid solver has succesfully been tested for advection of a square wave and evolution 
of a three dimensionnal Sedov blast wave. Finally, the MHD code gives good results on 
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various tests such as slow, fast and Alfven waves as well as an MHD shock-tube problem. 



To modify this TVD code to solve equations ([2]) through ([H]), we include an additional 
fluid variable and use exactly the same routine. The flux vector associated with the 
conservative variables becomes (for an update along 'x'), 

/ PVx \ 

fyvl+p + 0.5B'^/4:7r + Per - BI/Att 
pVxVy - B^By/A-K 
pv^v^ - B^BJA-K 
{etot +P + 0.55V47r + Per)v^ - B^B.v 



\ 



I 



with Per = (7cr— l)ecr- The freezing speed becomes c = max\\Vx\^\/ iciV + IcrPcr + -B^/47r)/p] 
and the timestep is reduced to dt = 0.8Ax/max[{\vx\, \vy\, \vz\) + ^/ {'jp + "jcrPcr + -B^/47r)/p]. 
In this way, we recover the original MHD TVD method when pcr = and otherwise take 
into account effects of cosmic rays using the same TVD routine. 

The remaining source term —p^r^.v is related to the pressure work during expansion 
or contraction and is easy to implement. Following the general philosophy of the code, we 
discretize it at second order accuracy using dimensional splitting for multidimensional runs 
and Runge-Kutta to reach second order temporal accuracy. It leads to: 



Fcri 



2Ax 



This contribution is finally included in the energy update at the beginning of each one- 
dimensional hydrodynamics step. In the following sections, we describe several tests of this 
2-fluids code that we have performed. 



3.1. Linear test: propagation of a sound wave in a composite of cosmic rays 

and thermal gas 

The first simple test is the propagation of sound wave in a medium with cosmic rays 
and plasma. The adiabatic wave speed is given by 



IPO + IcrPcrO ^^-^ 

Po 

with Po and Pcro the initial non-perturbed quantities. In order to trigger an eigenf unction, 
we need to satisfy the following relations between the field perturbations. 

So 6v , , 

- = 10 

Po Cs 



= 7cr-, (12) 
PcrO 

with 5v, 6p, 6p and 6pcr the perturbations. For our test, we take for the equihbrium quantities 
vq = 0, Po = 1, Po = 1 and Pcro = 1- We then perturbate the velocity with a sine of amphtude 
6v = 10'^ and wavelength 0.5. The other quantities are perturbated following [TOl [11] and 

Figj2] shows the results after a propagation during one period for 128 grid points. The 
result is in good agreement with the analytical solution and we obtain the same level of 
accuracy achieved in a pure hydrodynamical simulation (without cosmic rays). The slight 
smoothing of the extrema is due to the slope limiter which prevents the code from introducing 
spurious oscillations. 



3.2. Non Linear test: Riemann shock-tube problem for a composite of cosmic 

rays and thermal gas 



A more challenging t est is the Riemann shock-tube problem. The standard problem (jSod 



19781 ; iHawley et al.lll984j ) involves a polytropic gas starting with a state of high pressure and 
high density in the half-left space and a state of low density and low pressure in the half-right 
space. It leads to 5 regions with different fluid states separated by the head and tail of the 
rarefaction wave, the contact discontinuity and the shock. The interesting point is that one 
can derive the analytical solution using the Rankine-Hugoniot conditions. 

However, in our case the composite of cosmic rays and thermal gas is not a polytropic 



fluid and this solution doesn't apply. This problem has been solved by iPfrommer et al. 



(120061 ) and we use here their analytical solution. Our ID initial conditions are close to 
the ones used in their article with a left-hand state (L) and a right-hand state (R) in the 
simulation box. They are given (using an appropriate system of units) by, 

l<x<1.5 1.5<a;<2, (13) 

PL = 1 PR = 0.2, (14) 

vl = vn = 0, (15) 

PL = 6.7 X 10^ PR = 2.4 X 10^ (16) 

PcrL = 1.3 X 10^ PcrR = 2.4 X lOl (17) 

(18) 
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The sound speed is therefore Csl = 537 and Cs/? = 60. We run a simulation with 1024 grid 
points until t = 4.4 x 10~^ so that the shock front has propagated on an important fraction 
of the box length L = 1. 

Here again, there is a good agreement between the simulation results and the analytical 
prediction (see FigJSj). The transi tions between th e 5 states are well situated. The shock 



is resolved using few cells. As in iPen et al.l ( 120031 ) some variables have a slight overshoot 



in the first postchock cell but this doesn't affect the other subsequent cells. The contact 
discontinuity is slightly smoothed by the relaxation solver. To conclude, the accuracy is 
similar to the accuracy obtained in the 1-fiuid shock tube test and the implementation of 
cosmic rays is successful. 



4. Anisotropic transport: heat conduction and cosmic-ray diffusion 



In the presence of a magnetic field, the heat conduction in a plasma takes the form 
de 



dt 



V. Kubb.VT 



K^{I-bb)VT 



(19) 



with K± the perpendicula r conductivity an d the parallel conductivity and b the unit vector 
along the magnetic field (iBraginskiil 1 1 9651 ) . We will focus here on this equation, but one has 
to keep in mind that the diffusion of cosmic rays has a similar form 



dt 



V. {Dnbb.Vecr) + V 



D^{I-bb)Ve,r 



(20) 



with D_i the perpendicular diffusion coefficient and D\\ the parallel diffusion coefficient. This 
equation is therefore solved by the same subroutine. 

In cluster of galaxies the ion giroradius is much smaller than the mean free path between 
particle collisions and therefore the perpendicular part could be neglected since <^ k\\ 
(and Dj^ ^ D\\ for cosmic rays). The conduction is highly anisotropic, primarily along the 
magnetic fiel d, and mainly due to electrons. We adopt here the Spitzer value Hs for the 
conductivity (jSpitzeii 119621 ) of an ionised plasma with the Coulomb logarithm InA set to a 
typical value for clusters: 



Kll 



knT 



5 keV 



0.01 cm" 



Ks = 9.2 X lO'^VA^B 

For the parallel diffusion coefficient of the cosmic rays, we set 

Dn = 10^9 cm\s-\ 



37 
MA 



2 —1 

cm .s 



(21) 



(22) 
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4.1. Implementation and tests 

4.1.1. Centered asymmetric method: advantages and drawbacks 

The first metliod we implemented uses the so-called centered asymmetric differencing. 
It is the most natur a l conservati ve discretization and it has been shown to give good results 
by Parrish fc Stone (j2005 , 2007 ). The idea is to compute the heat flux F on each face and 



then to evolve the energy using an explicit time integration. We will consider only two 
dimensions but the generalisation to three dimensions is straightforward. The update of the 
energy is, 



„n+l n pn _ pn pn _ pn 



At Ax Ay 

This is the starting point for any conservative methods, now the problem is to evaluate 
the face-centered flux. The flux at time n and position {i + 1/2, j) is given by (see FigH]), 

(24) 
(25) 
(26) 

(27) 

(28) 

(29) 





dT - dT 
= b.Mb.^+by-), 




bx 


~ ^x,i+l/2,j^ 
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dT 


T T 




dx 


Ax 




by 


by,i,j-^i/2 + by,ij+i/2 + by^i+ij_i/2 


+ ^?/,i+lJ+l/2 


4 




dT 


^i+i,jr'+i ~ + Tij^i — Tij 


-1 


dy 


AAy 




component of the temperature gradient and the mag 


netic field are 



{i + 1/2, j), however the y components need to be extrapolated (overline). The time step is 
choosen to ensure linear stability 

/ Ax^ \ 

At = 0.45 min , (30) 

V2^cond/ 

^cond = (7-l)-^- (31) 



This m ethod is fast, efficient and accurate. However, as highlighted by lSharma &: Hammett 



(120071 ). this method is not positive definite. Indeed, it could lead to negative temperature 
in presence of large temperature gradient. An easy way to see the problem is to notice 
that in the flux expression Tj+i j+i and Tj appear with positive signs. So if one of this 
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temperature is a lot larger than all the others then nothing could balance this very large 
negative contribution and the energy e^^^ could become negative. This problem is due to 
the spatial discretization itself and not to the explicit scheme used for time integration. 
The transverse temperature gradient is not computed from the same origin as where the 



energy is taken, th is is the heart of the problem. An implicit scheme ( ISovinec et al.l 12005 



Balsara et al.l 120081 ) could therefore also suffer from the same negative temperature issues. 
One could indeed imagine the same situation as before but where the very large temperature 



stays relatively constant until the time step t""*"^. Then, the energy e"~j 



Ti+ij^i or Tij_i 
could also become negative with an implicit scheme. 



n+l 



There are two methods to make the scheme positive: the fir st one consists in limi t ing th e 
transverse gradient This idea is described in full detail in ISharma fc HammettI ( 20071 ). 
The second one consists in another discretization of the problem and is described in the 
following section. 



4.1.2. Positive anisotropic heat conduction: the flux-tube method 

We present here a new positive method for anisotropic conduction. It is based on 
a physically motivated discretization which treats anisotropic conduction as a ID diffusion 
process along the field lines. One could indeed simplify the discretization by considering only 
one thin magnetic flux-tube containing {i,j) and by calculating the temperature gradient 
and energy flux directly along this flux tube. Using V.-B = 0, the anisotropic conduction 
equation 



could be rewritten as 



I - B.v{^kVT), (33) 



b.vC^b.VT], (34) 



or, 



Bdt \B 

where b = B/B. One can see the apparition of the derivative along the magnetic field, and 
the term 1/5 to satisfy magnetic flux conservation. Because B.A = constant, where A is 
the cross-section area of a fiux-tube, the 1/B term can be thought of as representing the 
cross-sectional area A. 
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We now define s as the curvilinear abscissa along the field lines. The origin of this 
curvilinear abscissa is choosen to be s = at the grid point of interest. In order to 
compute derivatives, we consider variations over a length As = Aa; along the field line. The 
derivative of a function in becomes (/+as/2 — f-As/2)/^s and the same function 
derived in +As/2 gives (/+as — /o)/As. 



The discretization along the fiux-tube is therefore (see FigjS]), 



Bo 



At As 



B ds ) _^_^^/2 \B ds ^ _^^/2 



(35) 



= ' (36) 



9s J +AS/2 



^ (37) 



_ ^||,-As + ft^ll.O , . 

«||-As/2 - ^ , I3»j 

/t||,+A. + l^\\,0 , . 

f^\\,+As/2 = , (39) 

B_As + Bp 

B-As/2 - , (4Uj 

B+As + Bq 

B+As/2 - , (41) 

with the subscript indicating the curvilinear abscissa where the value is computed and the 

overline meaning that the value is not directly known and has therefore to be interpolated 
from grid point values. 

The next step is to estimate the position X[s) corresponding to s = ±As in the cartesian 
grid, which is done at second order, 

X{s) = bos + 0.5s\bo.V)bo, (42) 



b 



(6o.v)^o = + (43) 

(^o.V)6„o = b^J-^t^±l^^^^^ (44) 

The final and fundamental step is the interpolation of the temperature at the curvi- 
linear abscissa s = ±As. This will determine the accuracy of the method, as well as the 
positivity of the scheme. For this purpose, we decompose the 2D interpolation into a series 
of ID interpolations that are done with the second-order Lagrange interpolating formula. 
An important point, is that whatever the position we consider, we only interpolate using 



bx,i—l,j 7 bx^ij — l 
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and the 8 surrounding points. This is more convenient for the boundary conditions. 
Unfortunately, second-order interpolations are not guaranteed to stay in the range defined 
by the 2 extrema of the 9 considered point. Allowing such overshoot could create oscillations 
and negative temperature. We therefore saturate the interpolation to the extrema of the 9 
considered point, in order to allow positivity of the scheme. One drawback is that we loose 
accuracy near extrema, but this is unavoidable in order to get physical results. We will also 
see that the resulting amount of artificial diffusion is small. Finally, the norm of the magnetic 
field and the conductivity are interpolated only at first order for speed and therefore don't 
need to be saturated. One could interpolate at higher order for better accuracy. 

It is interesting to note that one could rewrite the update of the energy as, 

To"+i = + a f ^-^^/f-^^ + -^+A./2T+A. _ \ ^ ^^^^ 

A _ ^||,-As/2 Bp 

A-As/2 - — , (46) 

, _ ^||,+As/2 Bp 

A+As/2 - — -B , (47) 

f^lO -D+As/2 

" = ^cond,o^(^-A«/2 + ^+As/2), (48) 
Po Kb 

It means that the temperature Tq evolve by a fraction a toward the arithmetic average 
of T__As and T^as- We therefore choose the time step like in the precedent method that is to 
say. At = 0.45 min {Ax'^ /{2D^^j^^)). Using this time step and computing the norm of the 
magnetic field by |40] and SH it guarantees that a < 1 and prevents from any overshoot of the 
average. Since the interpolated temperature T_as and T^as are between the extrema of the 
neighbors of Tq, it means that no oscillations could appear! We have therefore implemented 
a positive fiux-tube scheme for anisotropic conduction and diffusion. 



4.1.3. Diffusion of a step function 

The first test we run is the passive diffusion of a ID Heavyside function. The goal here 
is to check if the code solves well the diffusion along straight magnetic field lines. In a second 
test we will check how well the code follows curved magnetic field line. We start here with 
the following conditions. 



p = 1 , 6 = 1 everywhere 



(50) 
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e = 1 for x< 0.5 (51) 
e = 2 for 0.5 < X < 0.75 (52) 
e = 1 for x> 0.5 (53) 

We then use a constant conduction coefficient, -D^qj^^^ = 1 so that the solution is ana- 
lytically tractable. Our 100 grid points simulation is ran until t = 2.8 x 10^^. For one step 
of size Ae and mean cq situated at the location xq the analytical solution gives, 

/ ,N Ae J X - xq \ , ^, 

e(x, t)=eo + —erf . (54) 

V V ^^cond^ y 

The comparison in Figj6] shows a very good agreement between the simulation and the 
analytical solution since we cannot differenciate them. 



4-1 .4- Anisotropic conduction in circular magnetic field lines 



A more c hallenging test iiivolves passive anisotropic diffusion along circular field lines, 
as proposed in lParrish fc Stond (120051). The idea is to consider an initial hot patch embedded 
in circular magnetic field lines. The heat should then diffuse along the field lines but not 
across the field. We start here with the following initial condition: 



p=l 
y-0.5 



for 
for 



r 

X - 0.5 



r 

10000 



for 
for 

1 otherwise, 



< X < 1 and < y < 1, 
< X < 1 and < y < 1, 

< X < 1 and < y < 1, 

0.7 < X < 0.8 and 0.49 <y < 0.51, 



(55) 
(56) 

(57) 
(58) 
(59) 



with r = ^{x-0.5)^ + {y-0.5y 



The 100 by 100 simulation is run with D^^^^ = 1 using our flux-tube method as well 
as the st andard centered asyrnmetri c method. We also run the Van-Leer-limited implemen- 
tation of ISharma fc HammettI (120071 ). We present in FigJ7|the temperature profiles at t = 0, 
t = 0.0225, t = 0.0675 and t = 0.18 for these three methods. 

In all three methods, the heat flux follows the circular field lines and tends toward a 
stationnary solution without any angular gradient of temperature. When there is no perpen- 
dicular conductivity, the analytical stationnary solution is obtained by energy conservation: 
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Gstat = 128.3 everywhere inside the shell. However, second order truncation errors add some 
artificial perpendicular diffusion. As a consequence the radial profile which was initially 
a Heavyside (as in the preceding test) diffuses. The consequence is that the maximum is 
lowered and the radial profile is smoothed. 

If one uses the standard asymmetric differencing, the dramatic consequence is that it 
leads to negative temperatures, even after a long run time. These negative temperatures are 
shown as a white inner and outer circle with dotted contours in FiglTJ left column. This is 
a very important problem. On the numerical point of view, it indicates that this scheme 
could overshoot the extrema and therefore create some spurious oscillations. On the physical 
point of view, it means that heat can flow from lower to higher temperature. Moreover, while 
coupling with the MHD solver, it could lead to negative temperatures, create an imaginary 
sound speed, and lead to unphysical results. 



All these points have been discussed in detail in ISharma fc HammettI (120071 ). They 
found that this problem arises in presence of strong temperature gradient perpendicular 
to the magnetic fleld. This is why they proposed a limited version of this asymmetric 
discretisation, in which they limit the perpendicular temperature gradient. However, as 
they already mentionned, the perpendicular diffusion becomes important if one uses limited 
methods. This is obvious, in FiglTJ middle column. For example, in the last line, one 
could see that the radial dispersion is larger than in the asymmetric method. Moreover, the 
maximum has been decreased by a factor of 1.5. 

On the contrary, our method combines two advantages of the two other methods. As 
presented in the right-hand column of FiglTJ the temperature always stays between the initial 
extrema but keep a low level of perpendicular diffusion. We are now going to estimate this 
perpendicular numerical diffusion using a test especially dedicated for this purpose. 



4.1.5. Accuracy of the method: Sovinec test 



In order to compare the accuracy of different methods, it is interesting to know what is 
the artiflcial perpendicular diffusivity of a scheme. Indeed , some applica t ions c ould require 
a large ratio of the parallel to perpendicular conductivity. ISovinec et al.l (120051 ) have devel- 
opped such a test. We will therefore run this test for our method and co mpare our perpendic 
ular a rtiflcial diffusion with the Van- Leer-limited method presented in ISharma fc Hammett 
( 120071 ) as well as the standard asymmetric method. 



The idea is to consider the full heat equation [19] in 2D with both a perpendicular and 
an anisotropic part. We also add in this energy equation a heating source term Q{x,y) = 
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Qq X cos{kx)cos{ky). The equation to solve becomes, 

K^{I-bb)VT +Qix,y). 



de 
di 



V.K||b6.VT + V. 



(60) 



The analytical stationnary solution could be computed in the case of pure isotropic conduc- 
tion = K±) with a constant conductivity. It is given by 



T{x,y) 



2K±k' 



-,cos{kx)cos{ky). 



(61) 



As in ISovinec et al.l (120051 ) . we consider a fixed magnetic field satisfying B.'VT = 0, so 
that the previous solution still apply. Taking into account artificial diffusion, the solution in 
the center becomes T(0, 0) = Qq/ [2A;^(k^ + ^num)]- The artificial diffusion could therefore 
easily be deduced from the central temperature. 



Our initial conditions for 0.5 < x < 0.5 and 0.5 < y < 0.5 are. 



p 


= 1 


k 


= vr 


Qo 


= 2n^ 




= 0, 


T{x,y) 


= cos{nx)cos{ny), 


Bx 


= cos{'iTx)sin{7iy), 


By 


= —cos{7cy)sin{TTX 



(62) 
(63) 
(64) 
(65) 
(66) 
(67) 
(68) 



We also choose Tt 



bound 



for the fixed boundary conditions. 



Unlike ISharma &: HammettI (120071 ) and ISovinec et al.l (120051 ) . our goal here, is to esti- 
mate the numerical diffusion in the case of a pure anisotropic conduction {k± = 0). This 
numerical diffusion increases with the parallel conductivity with k± = 0. We obtain Knum 
by running simulation to steady state and setting Unum = ^/T{0,0). Using this method, 
we determine the ratio Unum/'^w as a function of the resolution {dx/L = 0.01, dx/L = 0.02, 
dx/L = 0.05 and dx/L = 0.1) for the different implementations of the conduction. The 
results are presented on FigJHl 

The first point is that all the methods converge towards lower numerical diffusion with an 
order of convergence of ~ 2 (i.e., Unum/ dx^). The least diffusive method is of course the 
standard asymmetric method. This method reaches a ratio of Unumj = 10~^ for a resolution 
dx/L = 0.01. However, we have already noted in the precedent part that this method could 
lead to unphysical results and may therefore not be suitable for applications with large 
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tempe rature gradients (like in presence of shocks). The method of ISharma &: Hammett 



( 120071 ) circumvents this problem but the perpendicular diffusion is a factor of ~ 7 more 
important than in the standard method. On the contrary our new method is only a factor 
of ~ 2 more diffusive than the standard method but doesn't lead to negative temperature. 



4.2. Conclusion: comparison of the three methods for anisotropic conduction 

or diffusion 

We summarize in Table [T] the properties of the three diffe rent methods studied i n this 



section. One of the most important property emphasized by ISharma fc HammettI (120071 ) 
is to know if the solution remains bounded in the initial range of temperature. Indeed, 
this is essential to guarantee physical results and stability of the scheme. Unfortunately, 
the standard asymmetric method doesn't share this property. It still could be used in 
presence of smooth temperature field taking advantage of its speed (4 times faster than the 
MHD solver) and accuracy {i^num/ = 10"^) but has to be avoided in presence of strong 
temperature gradient and chaotic magnetic field. 

From our knowledge, only two methods for asymmetric conduction (or diffusion) are 
positive definite. In the first meth od, an asymetric discretiza tion is used but the transverse 



temperature gradients are limited (ISharma fc Hammettll2007l ). This method is almost as fast 



as the precedent one (three time faster that the MHD solver) but the limiter increases a lot 
the perpendicular diffusion (K„„rrt/'^|| = 7 x 10~^). The second method, from this article, is 
based on a physically motivated discretization along the magnetic flux tube. The accuracy 
then turns out to be better [KnumI ^^\\ = 2 x 10~^) but it is a little bit slower (although 
still one time and half faster than the MHD solver). In order to allow the reader to judge 
which problem size can be realistically treated we give here an indication of the cpu time 
and the memory consumption for a 3D run with 100^ grids on an AMD Operon 1.8 GHz. 
A double time step consisting in two calls to the transport subroutine and two calls to the 
MHD+cosmic rays routine takes about 30s and uses about 500 MB of memory. We are 
currently working on a parallel version in order to make larger runs. 

Since this method is a non-conservative method, we have also estimated the average 
fraction of energy lost per time step in the circular conduction test. These losses are limited 
to about 10~^, which is small considering that the magnetic field is strongly curved and the 
temperature falls by a factor of 10^ in few cells. Finally, it is worth noting that even higher 
level of anisotropy could be reached by implementing a higher order method. This could be 
easily done by interpolating the temperature field at higher order. 
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In the present and past section, we have shown that cosmic-ray tests and passive con- 
duction tests were succesfull. We now move to active conduction tests which involve couphng 
between the MHD solver, the cosmic-ray solver, as well as the anisotropic transport solver. 



5. Buoyancy instabilities 

In this section, we focus on buoyancy instabilities in a stratified atmosphere. We present 
here two applications of our code which serve both as a test of our two-fiuid code with 
anisotropic transport as well as a physical study of the cosmic ray magnetothermal instability 
(CRMTI). 



5.1. Physical background 



The cosmic ray magnetothermal instability (CRMTI) (jChandran &: Dennisll2006l:lDennis fc Chandran 



2007 ) is a buoyancy instability that is similar to the Parker instability (jParkei 



1966 



Shu 



1974 : iRyu et al.ll2003l ) since it involves magnetic fields and cosmic rays. However, in contrast 
to the Parker instability, the CRMTI involves anisotropic thermal conduction and allows for 
a temperature gradient in the equilibrium, both of which are relevant for understanding 
buoyancy instabilities in clu sters of galaxies. The CRMTI is very s i milar to the magne- 
tothermal instability (MTI) Jsalbuslboool . I2OO1I : [Parrish fc Stonelboosl . l2007h . except that it 
involves cosmic rays, and, m agnetic buoyancy if (3 = 87rp/B^ is n ot large (as in the analysis o f 



Dennis fc ChandranI (120071 )). The stability criterion is given by lDennis fc ChandranI (120071 ) 



nksdT/dz + dpcr/dz + dcmag/dz > 0, 



(69) 



with e 



mag 



B^/8ti the magnetic energy density. As discussed by IChandranI (120051 ) : IChandran fc Raseral 



(I2OO7I ). the CRMTI may lead to convection in galaxy cluster cores, since central AGN pro- 
duce jets and centrally concentrated cosmic rays. Such convection may play an important 
role in transferring AGN power to the intracluster medium and helping to solve the "cooling 
fiow problem." 

To understand the physics of this instability, let's take a simple example where this 
instability applies. Consider a magnetized plasma plus cosmic-ray stratified atmosphere 
initially at equilibrium with a negligible gradient of magnetic field strength and temperature 
but a negative gradient of cosmic-ray pressure {dpcr/dz < 0). Consider also a vertical gravity 
along "z" g = —ge^ and a horizontal magnetic field along "y", Bq = BoSy. If one perturbs 
this medium by pushing upward a fiuid parcel {6vz > 0, with 6 the difference between the 
perturbed state in the bubble and the equilibrium state in the surrounding), the magnetic 
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field lines will be distorted [SB^ 7^ 0). As a consequence of the anisotropic transport, the 
cosmic ray pressure and plasma temperature will be smoothed along the perturbed field 
lines, that is to say the temperature and cosmic-ray pressure of the upwardly displaced fluid 
parcel will be the same as in the initial equilibrium. Therefore, the cosmic ray pressure in 
the upwardly displaced fluid parcel will be larger than the one of the surrounding medium 
{Spcr > 0) and the temperature will be almost the same {ST f» 0). Requiring total pressure 
equilibration with the surrounding medium, it implies that the parcel density is lower than 
the one of the surrounding medium {6p > 0) to compensate for the high cosmic ray pressure. 
As a consequence, the parcel moves upward faster, and the medium is convectively unstable. 
This instability could even be amplified if the temperature gradient is negative since in this 
case we would have 6T > 0. The latter occurs in the magnetothermal instability (MTI), 
where dp^r/dz = but dT/dz < Jsalbuslboool ). 



An important assumption in the above analysis is that Bzo = 0. If instead Bzo 7^ 0, 
then there i s an equilibriurn heat flux in the z direct ion, which can further contribute to 
instabilities (iQuataert II2OO8I : iParrish fc Quataertll2008l ). However, such "heat-flux buoyancy 
instabilities" are not considered further in this paper. 



The linear analysis fo r the CRMTI has been done by IChandran &: DennisI ( 120061 ) and 
Dennis fc ChandranI (120071 ). The dispersion relation is given by 



+ 00' 



dhipo 

dz 



(70) 



ekyA2u' + vl) - [kl + kl) { / + in' + v\)g^] 



dz I 



+ klvl 



72722 2 , n 2 , 7 2\ I 2, 2 dlnpo 

-k kyV^u +{k^ + ky) [g +u 



with uj the frequency, k 
Alfven speed and 



.,ky,kz) the wave vector, v\ = Bq/^Ahpq) the square of the 



u 



PO + ir] ^ PcrO Icr^ 



V 

V 



Po u + ii] ' po LU + ir] 

k'^D 
y cond' 



(71) 

(72) 
(73) 



The eigenfunctions are given by, 

5p = 
SB = 



iSvz dpo ^ k.Svpo 



LU dz 



LU 



iSvz dBf) kyBoSv k.SvBg 

I • 



LU 



dz 



LU 



LU 



(74) 
(75) 
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5p 



5pc 



6vr 



ijj dz 
i5vz dpcro 



H.F 



dz 
C.E 



A.E - 
A.F-C.H 

- A.E 



+ 



5Vy 



UJ + 17] UJ 
k.dv-fcrPcrO 



(76) 
(77) 
(78) 
(79) 



with A 
F 



-l\LO 



kyvl - kliu^ + vl)], C = -Kg + iKk,{u'^ + v\), E = -i{uj^ 

i k^kyV?. For a d etaile d discussion of this instabihty see 



kyQ and H 



Chandran &: Dennis! (120061 ): iDennis fc ChandranI (120071 ) . Let's now study this instabihty 



using our new code. 



5.2. An interesting limit: the magnetothermal instability (MTI) 



The easiest limit of this system is the adiabatic convective instability (no cosmic ray, no 
magnetic field, no conduction) which obeys the Schwarzchild stability criterion {dS/dz > 0). 
As a preliminar test we have run a simulation and compared the evolution of the eigenfunction 
with the analytical one. The agreement between simulation and linear analysis is good. We 
do not present here the results since it has been abundantly studied in the literature. 

A more complex limit arises in dilute plasma when one takes into accou nt anisotropic 
conduction. Indeed, in this case, the stability criterion becomes dT/dz > (IBalbusI l2000l ) 
and the medium could therefore be unsta ble even if the eri t ropy gradient is positive. Since 
this instability has already been studied in jParrish fc Stond (120051 ). our main purpose here is 
to test our code by comparing the numerical s olution with the li near analysis. We therefore 
use a very similar test case to that studied by jParrish fc Stond (120051 ). The main difference 
is that we trigger only two modes so that we are able to compute the exact linear solution 
(eigenfunctions) . 

Our initial conditions for the magnetothermal instability test are therefore a vertical 
eq uilibrium state for — .05 < y < 0.05 and —0.05 < z < 0.05 (we use here the same axis as 



m 



Chandran &: Dennis! (12006! ) even if it's a 2D run): 



p{z) 
9 



lip 



m 

(81) 
(82) 
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B 



with pq = pq = Tq = 1 ui appropriate u nits and 
By choosing Hp = 1.5 and Hp = 1, as in 



B„ (83) 

(84) 

2 X 10^ (t o be in the high beta hmit). 



Parrish &: Stond (120051 ) . we ensure that the entropy 



gradient stays positive but the temperature gradient is negative. The conduction coefficient 
is -DpQj^(-^ = 6.8 X 10^^. The box length L = 0.1 satisfies the condition L < 0.1min{Hp, Hp) 
which is essential to be able to apply local analysis. 

The boundary conditions need to conserve energy and also to not break artificially the 
initial equilibrium. We therefore choose periodic horizontal boundary conditions and reflec- 
tive vertical boundary conditions. The implementation of reflective boundary conditions in 
presence of gravity turns out to be non trivial since the last cells are not in equilibrium if 
one uses standard reflective boundary conditions. In order to deal with this problem, we 
assign values to the ghost cells by assuming reflectional symmetry for scalars and reflectional 
antisymmetry for vectors, including gravity. We also remove the artiflcial diffusion of the 
last active cells that is normally added explicitly as part of the TVD method. In this way 
we could fuUflll the two requirements: energy conservation and equilibrium. 

We perturb the initial equilibrium with the superposition of two eigenf unctions with 
same ky = 2'k/L but opposite {k^ = +27c/L and k^ = —2tt/L) so that the vertical velocity 
cancels along the horizontal boundaries. We take the amplitude of the velocity perturbation 
of each mode to be 6vz/cs = 1.55 x 10^^ so that all the perturbed flelds stay in the linear 
regime. Taking such a small value is very important since we will see that the relative 
amplitude of the magnetic fleld fluctuations are a factor ^ Csk/a ^ 200 times larger than 
the one of the velocity. To compute the eigenfuncti on we take the limit Pcr = and = 
of the eigenfunction of IChandran fc Dennis! (120061 ) described in the last subsection. We 
flnally study the most unstable mode which is convective and exponentially growing. For 
this purpose, we run a 2D simulation with 200 by 200 grid points. The expected growth 
rate is a = 0.22. Note that to obtain the same results as Parrish fc Stone (12005) one has to 
trigger only one eigenfunction with k^ = but, in this case, boundary conditions then excite 
a lot of different and uncontrolled modes. 

The results in Figl9]show good agreement with the analytical result. We however note 
some slight deformations mainly due to the slope limiter as already mentioned in earlier 
sections. For example, one could clearly see the clipping of the maxima of 6vz- We have 
checked that without slope lirniters t hese deformations do not appear. We therefore recover 
the results of iParrish fc Stond (120051 ). concerning the growth rate. We have presented a test 
for active anisotropic conduction which turns out to be very strict and give no mercy to any 
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approximations in any part of the scheme. We have noticed that a simple precision run or 
a too violent slope limiter quickly destroy the shape of the sines. 



5.3. Solution of the full system: the cosmic-ray magnetothermal instability 

(CRMTI) 

This code is especially dedicated to the study of this instability. This part serves as 
a cros s-validation between our code and the analytical predictions of IChandran &: Dennis 
(120061 ) . We use here the same kind of initial condition as in the study of the magnetothermal 
instability. The difference is that we have now cosmic rays and we therefore take, 

z . 



Vc 



Pcro(l 



(85) 
(86) 



Instea d of choosing the same density and pressure gradients for the plasma as lParrish fc Stone 



(120051), we prefer to study a different situation where the gradient of entropy and the gradi- 



ent of temperatu re are positive. In this case, both the Schwarzchild stability criterion and 



the iBalbud ( 120001 ) stability criterion are satisfied. However, we choose a negative gradient of 



situation is expected in the center of cluster of galaxies ( 


Chandran 


2005; 


Chandran & Dennis 


2006; 


Dennis & Chandran 


2007; 


Chandran & Rasera 


2007; 


Rasera & Chandran 


2008). Since 



we are mainly interested in the role of cosmic rays, we take a n ul gradient of magnetic fiel d 
so that the magnetic field doesn't modify the stability criterion (IDennis fc ChandranI 120071 ). 



^Our initial conditions are inspired by the Perseus cluster at 50 kpc (IChandran fc Dennis 

20061 ). Namely, we take Tq = 4 keV, an electronic density ngo = 0.02 cm~^, an ion density 
of rijo = 0.9ne, a mean molecular weigh of /i = 0.6, a cosmic ray pressure pcro = Po and 
a magnetic field of -Bq = 1 /uC For these values, the conduction coefficient is ^q^q^A ~ 
2.5 X 10^° cm^/s and we choose the diffusion coefficient to be Dy = 10^^ cm? .s~^ (see 
Sectj2]). Concerning the gradient, we take a negative density gradient with Hp = 50 kpc, 
a positive temperature gradient with Ht = 200 kpc, and a negative cosmic ray gradient 
with Her = 50 kpc. We finally use 100 by 100 grid points for a simulation box length of 
L = 10 kpc. The amplitude of our two perturbations is Svz/cs = 1.44 x 10~^ and their 
wavelength is 10 kpc. We use the same boundary conditions as before and the expected 
growth rate is a = 9.5 x lO^'^Myr^-'^. 



Here again the results show good agreement with the linear theory (Fig lTOl) . The ex- 
trema of velocity are affected by the limiter. Little numerical errors appear on the temper- 
ature fluctuations because they have the smallest relative amplitude and are therefore the 



-24- 



most sensitive to numerical approximations (such as the one induced by the slope limiter). 
The cosmic-ray pressure is computed with the same accuracy as the plasma pressure. This 
validates our overall implementation. 

On the physical point of view, this simulation illustrates the scenario presented in part 
15.11 The phase of 5Bz is shifted from 'k/2 compared to the phase of 5vz-, since the field 
lines are distorted. The fluctuations of temperature 5T are very small. The perturbation 
5pcr is roughly in phase with 5vz and is amplified, which in turn amplifies 5vz- Even though 
the entropy and temperature gradients are positive (and the magnetic field energy gradient 
is nul), the convective instability is growing on a timescale of 10® yr which is of order of 
the cooling time near the center of cooling-fiow clusters. This suggests that cosmic rays and 
anisotropic transport might play an important role in cluster of galaxies. 



6. Conclusion 

This article could be viewed as a test guide for those who want to implement cosmic-ray 
and anistropic transport routines, which are essential ingredients to simulate cooling-fiow 
clusters. We indeed used many linear and non-linear tests for each physical ingredient as 
well as a new linear test for the full system. Our contribution could be divided into three 
parts. 



First, we showed that the TVD method of iPen et al.l (120031 ) can be used to evolve the 
cosmic rays and the plasma simultaneously. The shock tube problem for a composite of 
plasma and cosmic rays is an example of the successful implementation. Second, we insisted 
on the importance of having a positive implementation of the anisotropic conduction in order 
to ensure physical results even in the presence of sharp gradients. We therefore presented a 
new fiux-tube method which has two important properties: positivity and accuracy. This is 
important for clusters of galaxies since the conduction is highly anisotropic as opposed to a 
very diffusive scheme. Moreover, the random magnetic fields and potential large temperature 
and cosmic-ray-pressure gradients in cluster cores may cause negative temperatures in a 
non-positive scheme. Third, the linear regime of the cosmic ray magnetothermal instability 
(CRMTI) provides a new, sensitive test, of the overall implementation. The main interest 
is that this instability has a broader range of applications since the criterion is nkBdT/dz + 
dpcr/dz + dcmag/dz > 0. One interesting future application of this code concerns the cores 
of clusters of galaxies. Indeed, in these regions, negative radial gradients of cosmic-ray 
pressure may trigger convection. Moreover, possible large gradients of cosmic-ray pressure 
near the edges of X-ray cavities (cosmic-ray bubbles) require positive implementation of 
the anisotropic transport. It is worth noting that even though our discussion has focused 
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mainly on clusters of galaxies, the flux-tube method that we have developed could be used 
to simulate a variety of physical systems, including the interstellar medium, fusion plasmas, 
and supernovae-remnants. 

We are grateful to U.-L. Pen, T. Dennis, P. Sharma, W. Hammett, R. Teyssier, I. Parrish 
and J. Stone for their valuable comments and suggestions. 
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Method 


positive 


diffusion 


speed2D 


speed3D 


losses 


standard 


NO 


10-^ 


4.2 


4.6 


10-10 


limited 


YES 


7 X 10-^ 


3.3 


3.8 


10-10 


flux tube 


YES 


2 X 10-^ 


1.6 


1.0 


10-^ 



Table 1: Summary of the three second order methods for anisotropic conduction or diffusion. 



Note. — Columns indicate respectively: the method (standard asymmetric discretization, Van Leer 
limited or flux tube method), the positivity of the solution, the value of k„„„i/k|| in the Sovinec test, the 
speed relative to the MHD solver in 2D and 3D CRMTI instability test (see Sect lS.Sp and the fraction of 
energy lost per time step at the end of the circular conduction test. 
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Fig. 1. — Staggered grid used by the MHD solver. Density, momentum and energy are 
defined at the cell center whereas magnetic fields are defined on the faces. Fluxes are also 
computed on the faces and depend on the neighbours (squares). 



-32- 



0.0005 - 



o 

c 



k:) 



0.0000 



-0.0005 - 




0.0 



0.2 



0.4 0.6 

X 



0.8 



1.0 



Fig. 2. — Sound wave in a composite of cosmic rays and thermal gas propagated for one wave 
period. Plus signs represent the simulation cosmic-ray internal energy for 128 grid points 
whereas the continuous line is the analytical solution. 
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Fig. 3. — Shock tube problem for a composite of cosmic rays and thermal gas. Plus signs 
show results for a 1024 grid points simulation whereas the continuous hue is the analytical 
solution. The dotted line is the initial condition. The first graph shows the density profile, 
the second one the velocity profile, the third one the gas internal energy profile and finally 
the fourth one shows the cosmic-ray internal energy. 
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Fig. 4. — In the centered asymmetric method, the heat fluxes on each faces (bold arrows) are 
computed by projecting temperature gradients on the magnetic field. Temperature gradients 
computation require the knowledge of the neighbours (squares). 



Fig. 5. — In the flux-tube method, one considers a thin magnetic flux tube around Tq. The 
evolution of this temperature depends on the heat flux along the magnetic field which is 
evaluated on the two faces at a distance As/2 from the center: p^^^l"^^ These fluxes are 
deduced from the temperatures T^^* which are themselves computed using all the neigbhours 
(squares) . 
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Fig. 7. — Anisotropic diffusion of a hot patch (first graph from the top) along circular 
magnetic field in a simulation with 100^ grid points. Columns correspond respectively to, 
the standard centered asymmetric discretization, the Van-Leer limited methods and our new 
flux-tube method. Lines correspond respectively tot = 0.0225, t = 0.0675 and t = 0.18. The 
circles are the field lines confining the heat. Note the negative temperature (dotted contours) 
in the first method and the considerable perpendicular diffusion in the second method. 
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Fig. 8. — Measure of the perpendicular numerical diffusion Knum/^ii in the ISovinec et al. 



( 120051 ) test as a function of the resolution. The dashed curve is the result for the standard 
asymmetric discretiz a tion. The dash-dotted curve is for the Van-Leer-limited method from 
Sharma fc HammettI (120071 ). Finally, the continuous curve is for our new flux-tube method. 
For reference, the dotted line has a slope of 2. 
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Fig. 9. — Evolution of the eigenfunction in the code (plus signs) compared to the linear 
theory (solid line). The initial condition is the dotted line. The 6 graphs represent a slice in 
z = 0.0125 for respectively Svy/cg, Sv^/cg, SBy/Bo, 6Bz/Bq, 6T/Tq and Sp/po . 
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Fig. 10. — Evolution of the eigenfunction in the code (plus signs) compared to the linear 
theory (solid line, in most case covered by the plus signs). The initial condition is the dotted 
line. The 6 graphs represent a slice in 2; = 6.2 kpc for respectively, 6vy/cs, 6vz/cs, 6By/BQ, 
SB^/Bq, 6T/T0, 6p/po, Sp/po and Spcr/Pcro- 



